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Abstract 

An algorithm for generating optimal nonuniform grids for solving the two- 
body Schrodinger equation is developed and implemented. The shape of the 
grid is optimized to accurately reproduce the low-energy part of the spec- 
trum of the Schrodinger operator. Grids constructed this way are applicable 
to more complex few-body systems where the number of grid points is a 
critical limitation to numerical accuracy. The utility of the grid generation 
for improving few-body calculations is illustrated through an application to 
bound states of He trimers. 

Keywords: Quantum few-body calculations, Nonuniform grids, 
Optimization, Grading function, Faddeev equations 

Introduction 

The dynamics of few-body systems remains a robust field of research with 
many practical applications. A number of theoretical advances, coupled with 
increased computational resources, have lead to significant advances in both 
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the understanding of few-body processes and in the number of physical sys- 
tems that can be successfully treated with existing, well-tested methodolo- 



gies (see 



7 



9] and references therein). We are, for example, 
currently developing public source code and a graphical user interface for sci- 
entists interested in solving Faddeev equations numerically for a wide array 
of potential three-body applications. Challenging computations that have 
been accessible to only a small group of specialists will soon become elemen- 
tary and well characterized tools used by a large number of practitioners in 
different fields. 

Solving the two-body Schrodinger equation numerically is an elementary 
exercise in the case of smooth central potentials. Typically, the power of 
modern computers makes it possible to use even the simplest numerical ap- 
proaches to perform quite accurate calculations of the low-energy part of the 
two-body spectrum. Three-or-more-body calculations, however, usually re- 
quire greater attention to the details of numerical technique, as the required 
computer resources usually scale geometrically with the growth of dimension- 
ality of the problem, and optimizing any aspect of the solution representation 
leads to substantial computational savings. 

When solving bound state or scattering problems for few-body systems 
it is important to treat the states of two-body subsystems carefully, as such 
states represent the asymptotic boundary conditions for the corresponding 
few-body states. It is also important to minimize the computational cost of 
reproducing every asymptote of the few-body calculations. The key feature 
of the Faddeev approach to few-body problem is the asymptotic factorability 
of the solutions, so that grids constructed for the efficient solution of the 
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two-body problem can be immediately employed with considerable numerical 



advantage. In previous calculations 
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13| this optimization has been 



performed manually. The procedure, however, is very time consuming and 
difficult for an inexperienced user or a student. Therefore, we needed an 
automatic procedure of constructing an effective grid representation of the 
two-body subsystems. 

Methods of refining grids automatically are often used in solving vari- 
ous nonlinear evolution equations (hydrodynamical equations, for example) 
to reproduce discontinuities and other singular features of the solutions. In 
contrast, our goal is to create a software package specifically designed to solve 
the quantum-mechanical few-body problem, fully exploiting the features in- 
trinsic to the physical problem and the numerical techniques to achieve high- 
performance of the resulting code. We therefore needed a solution which is 
on the one hand more specific to our problem, and on the other hand allows 
us to construct the grids on the base of some clearly understood physical 
and mathematical principles, with particular emphasis on reproducing the 
low-energy part of the two-body Schrodinger operator for subsequent use in 
more complex few-body calculations. 

In this paper we describe and implement a practical nonuniform grid 
suitable for reproducing the low-energy part of the Schrodinger operator for 
two-body systems. When applied to systems of more than two particles, this 
grid permits a several-fold reduction in the number of grid points required 
for a desired level of numerical accuracy. Equally important, the procedure 
of constructing the grid is automatic and requires only minimal interference 
from a user. 
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1. Basic principles 

We are solving the Schrodinger equation for a two-body system with cen- 
tral potentials using S 3t2 or S 5:3 Hermite splines and collocations at Gaussian 
points. We are constructing a nonuniform grid from the requirement of uni- 
formity of the numerical error over the entire range where the solution is 
constructed. 

Let us start from the radial Schrodinger equation in atomic units (a.u.) 
chosen so that H = e = m e = 1 

where /z is the reduced mass of the two-body system, subject to the Dirichlet 
boundary conditions 

<p(0) = <p{x max ) = . (2) 

The right boundary of the interval x max is assumed to be chosen so that its 
influence can be neglected. Let Ajv, x = {0, x±, . . . , xn} be a grid constructed 
over the interval [0,Xn = x max ] consisting of N intervals, x is the map 
used to construct the nonuniform grid from a uniform grid. Particularly, 
Xi = x max x{i/N). The map x '■ [0, 1] — > [0, 1] is a smooth function growing 
monotonically over the interval [0, 1]. 

The criterion we choose for constructing x is the uniformity of numerical 
error over the whole interval [0, x max \. An approximate solution constructed 
with iV points deviates from the exact solution by 

(p( N \x) = ip(x) + residue at(x) 

At each sub-interval [a^Xj+i] the residue norm can be estimated as 

||residueAr(x)|| < Ci(x)\x i+1 - Xi\ k+1 + o(\x i+1 - Xi\ k+1 ) 
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where x G and k is the order of the spline [10|. The factors Ci(x) 

are determined by the properties of the equation and the chosen numerical 
scheme. We can optimize the distribution of the grid points so that the error 
| \Ci(x) \x i+ i — Xi\ k+l \ \ has the same order of magnitude throughout the whole 
interval. It is useful to treat the length hi = \xi — of the i-th interval 
as a continuous function h(u). The function h(u) is naturally related to the 
derivative of the map x( u ) 

hi = (xi-Xi^) = x max (x(^) - xC-j^-)) = ^j^x' C ^ //2 ) +Q (]^) ( 3 ) 



and 



So, for sufficiently dense grids the z-th grid step is determined by the slope 
of the mapping function. 

If e(t) is some non- negative function, termed the "grading function", that 
characterizes the local approximation error, then the integral 

= JT~f)d* (4) 

provides a monotonic map connecting the non-uniform grid with a uniform 
one. xi u ) is referred to herein as the "mapping function". If the grading 
function is chosen so that it peaks in the regions that are the most difficult 
to reproduce numerically, the inverse of the mapping function x l { u ) wm 
grow the most rapidly in the corresponding regions. The mapping function, 
therefore, will have smaller derivative. It will make the nonuniform grid 
denser where the function is more difficult to reproduce. Characterizing the 
difficulty of representing the object function, the grading function can be 
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linked to some derivative of the function being interpolated. For example, if 
using a linear interpolant, the error of approximation will be of the order of 
the second derivative of the interpolated function. A quantitative characteri- 
zation of the error and optimal properties of the grading functions have been 
studied in Ref. We shall use the approach of Ref. jl| for constructing an 
optimal map for solving Eq. [TJ 

Suppose we have a sample function <p(x) for which we seek a piece- wise 



polynomial approximant (p(x). It is shown in Ref. [1[ that in order to mini- 
mize the L2-norm of the error function || tp — (p ||^ 2 the grading function e(t) 
should be chosen as 

d k+1 2 

where k is the order of the polynomial approximant. According to [1] this 
grading function provides an asymptotically L 2 -optimal grid for a sufficiently 
large N, when the 0(N~ 2 ) term in the equation ([3]) can be neglected. 

There are a few complications with implementing this approach directly. 
The exact sample function <p(x) is unknown. Moreover, our goal is an ac- 
curate approximation of an invariant subspace of the two-body Hamiltonian 
which is characterized by several functions that are the eigenvectors of the 
Hamiltonian. Those complications, however, are not critical. 

If we redefine the norm to be optimized for a vector function composed 
of a set of eigenf unctions of the Hamiltonian {{pi(x), ip 2 (x), . . . , ip n (x)}, the 
optimal grading function takes the form 

(EI^T^)I 2 )^ • (5) 

m=l ux 



ex 



Although the explicit form of the functions ip m {x) is generally unknown, for 
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our purpose we can use numerical approximation of these functions obtained 
on non-optimal grids. The resulting suboptimal estimates can be refined 
iteratively. 

In the following sections we shall provide a few examples of how this 
approach can be realized in practice. 



2. Implementation 

Before describing the grid construction algorithm itself, we shall briefly 
outline the discretization procedure for the simplest case of zero angular mo- 
mentum I = 0. For discretization we use the orthogonal collocation scheme 
10 1 with the splines We shall seek for a solution of Eq. (JTJ an 



expansion in terms of B-spline basis in the spline space (A) constructed 
for a given mesh A 

M 
3=1 

Functions Bj(x) are constructed to satisfy the appropriate boundary condi- 
tions. Requiring the equation ([T]) to be satisfied exactly in the given set of 
collocation points Xk, k — 1, 2, . . . , M we reduce the equation to the (gener- 
alized) eigenvalue problem 

(-D a + V - E m B)f m = . 



Here D2, V and B are square matrices 



1 „« 



[Vh = (Vfr) + l ^-)B j (x i ) 



[B] tj = Bjixi) . 

The number of collocation points M - and the number of elements of the 
basis Bj(x) - depend on the chosen spline space and the number of grid 
points N. For £3^ splines and Dirichlet boundary conditions M = 2N. In 
the case of £5^ splines M = 3N + 1. 

Following the strategy we described in the previous section, first we dis- 
cretize the two-body Hamiltonian using some non-optimal grid Aq. A cubic 
mapping Xo( u ) — u3 1S , usually, a good starting point. 

Next, we construct a sequence of grids A,; = A Xi)N and the corresponding 
approximations to the lowest n eigenfunction (f m ^(x). Each subsequent grid 
A i+ i is constructed using equations (@| and (j5J) with approximate functions 
ip m (x) = <f m ^(x) obtained at the previous step. We repeat the process 
several times until the grading function is stabilized. 

Straightforward implementation of this approach, however, may look un- 
realistic. Indeed, we seek a solution in terms of piecewise polynomial func- 
tions of 5-th order. Evaluation of the grading function, however, requires 
evaluation of the 6-th order derivative of the solution. We, therefore, can not 
differentiate the spline expansion of the solution directly. Instead, we ap- 
ply the discrete analog of the second derivative operator B~ 1 D 2 three times. 
Each application projects the corresponding derivative function back into the 
spline basis and this way an approximate 6-th derivative of the solution can 
be obtained. 

In Fig. [1] we show an example of the sequence of the grading functions 
and the corresponding mappings Xi( u ) f° r t ne LM2M2 (He-He) potential jjjj 
and Ne-Ne potential (av5z+(3321) fit) [15j for S 3>2 splines (for I = angular 
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momentum). In both cases we have chosen the number of states n in Eq. <^ 
as the number of bound states supported by the potential plus one, so that 
the grid is optimized to reproduce all the bound states plus the continuum 
state closest to the threshold that satisfies the boundary conditions (T5]). (The 
He-He potential supports one s-wave bound state and the Ne-Ne potential 
supports three.) 

It takes only a few iterations for the grading function to converge. The 
grading function reaches its maximum in the potential well, rapidly falls 
at the repulsive core and slowly decays at larger distances. The inverse of 
the mapping function x _1 (w), according to equation (TjJ, rapidly grows at 
the values of u that correspond to the potential well, so that the mapping 
function xi u ) grows very slowly in a wide range of the parameter u. The 
grid steps, therefore, become much smaller in the important potential well 
region (see Eq. (J3])). 

In the next section we shall compare grids generated by the optimal map- 
pings with other empirical choices. 

3. Convergence of two-body bound states 

How practical is the construction? Can we expect any computational sav- 
ings in few-body calculations when using the optimal mappings compared to 
widely used simple power or exponential mappings? To answer this question 
let us study the convergence of the bound state energies for the two-body 
problem with respect to the number of grid points. 

To study the convergence properties of the numerical procedure it is natu- 
ral to represent the property of interest - the bound state energy for example 
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a) c) 




Figure 1: Grading functions and the corresponding mappings for He- He (left, a) and 
b)) and Ne-Ne potentials (right, c) and d)). The grading functions are maximal in the 
integration region ~ 1 — 10 a.u., as a result, the corresponding mapping functions have 
smaller derivatives ensuring smaller grid steps in the region of interaction. 
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- as a sum of two terms: the exact value and the grid- dependent numerical 
error 

E w = E + E err (A NjX ) . 

The numerical error depends on the number of grid points N, their distribu- 
tion defined by the mapping x an d the boundary conditions. Here we shall 
assume that the boundary conditions are chosen so that they essentially do 
not contribute to the error, or can be taken into account analytically. We 
shall consider the following grid shapes 

"Y n = H n 

A.power ) 

n _ e nu -l 
Xexp e n — 1 

for comparison with \opt defined by Eq. (jlj) and (jSJ) . To make this comparison 
we want to introduce a quantity which describes the character of convergence 
qualitatively. 

The collocation methods used in this work are expected to converge with 
the rate 0(iV~( fc+1 )) where k is the order of the spline. It is, therefore, natural 
to parametrize the numerical error by the inverse number of grid points -4. 
The error term E err (j^) does not behave regularly with the number of grid 
points and it is difficult to estimate it exactly. Its absolute value, however, 
should behave as N l +1 . To characterize the convergence obtained with a 
particular mapping \ we shall fit the numerical values obtained with the 
given number of points E^ as 

E^^E + Cj^. (6) 

The coefficient C characterizes the speed of convergence and E gives a more 
accurate estimate for the energy of the bound state. Both the speed of 
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convergence C and the extrapolated energy estimate E discussed below are 
obtained by least square fits of the expression (jSJ) to data sets similar to 
the one shown in Fig. [2j When fitting we gave bigger weights to the values 
obtained with denser grids. Having a finite number of points in each data 
set, we can extract estimates for E and C only with finite accuracy. As we 
have mentioned above, the error term does not behave monotonically, which 
makes it difficult to obtain very accurate estimates of E and C. E For the 
extrapolated energy estimate E we make sure that the error of the estimated 
value is comparable with systematic errors of different origins, such as the 
error introduced by an approximate boundary condition at the right end of 
the interval. 

To illustrate the meaning of the coefficient C, in Fig. [5] we show the 
convergence plots for the bound state energy of He2 calculated with three 
different mappings: Xo P t, xl xp an d Xpower- I* 1 Fig Ek) the fit of the bound 
state energy convergence for £3,2 splines is shown. The method of orthogonal 
collocations that we employ is not variational, and the error for the Hamilto- 
nian eigenvalues is expected to have the same order of accuracy as the wave 
function. The computational error, therefore, should scale as 0(N~ A ). We 



1 We see two sources of this non-monotonic behavior. First, the collocation scheme that 
we employ does not give an optimal variational estimate for the energy, therefore, when 
the grid is sparse, the deviation of the energy estimate from an optimal value is comparable 
with the error of approximation of the wave function. There is also a second source of non- 
monotonic behavior which influences even variational estimates of the energy and becomes 
evident for very dense grids: when increasing the number of grid points the spline space 
itself does not approach the invariant subspace of the Hamiltonian monotonically. This 
second effect, however, is much smaller than the first one. 



12 




a) b) 
Figure 2: Convergence of the He2 bound state energy for LM2M2 potential 

do observe this kind of scaling in Fig. (2^-b. The speed of convergence coef- 
ficient C corresponds to the slope of the lines in Fig. |2j It is important to 
note that as the coefficient C is calculated from the fits, its typical relative 
error is of the order of 50%, its estimate may vary depending on a particular 
grid sampling used in the fitting procedure. We thus will report only one 
significant figure in the estimates of the coefficient C, and deviations of the 
coefficient estimate within the same order of magnitude can be considered 
insignificant. Accordingly, C gives the desired qualitative measure of a grid 
quality. 

In Table [TJ we report the fitted values of the speed of convergence coeffi- 
cient. It is evident, that the error estimate for the optimal mapping is at least 
an order of magnitude smaller than for other considered grid choices. The 
result is not surprising, as can be seen from Fig. |3j Indeed, as the optimal 
mapping is not linear in logarithmic scale, simple power grids x™ OW er are n °t 
expected to produce a good approximation of the solution. The exponential 
mappings, on the other hand, can fit the optimal mapping more closely at 
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n ~ 8, and the estimates of the speed of convergence for n = 8 are minimal 
within the set of exponential grids. 

All the grids we demonstrated have been constructed for two-body states 
with zero orbital angular momentum / = (EqJT]). When constructing grids 
for actual use in few-body calculations the bound states with higher angular 
momenta should also be reproduced accurately. The centrifugal barrier for 
these states is usually a small perturbation compared to the two-body force, 
so that at least for a few lowest angular momentum states explicit inclusion of 
the rotationally excited states into the grading function may be unnecessary. 
Let us compare the spectra of rotationally excited Ne2 dimer obtained with 
the grid optimized for the s-wave with the spectra calculated with grids 
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Tabic 1: Estimates of the convergence speed C for the Hei and Nei s-wave bound state 
energies for different grid shapes 15 
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Table 2: Ne2 speed of convergence of the rotationally excited states. The optimal mapping 
Xopt has been constructed for the states with angular momentum quantum number l op t, 
the corresponding speed of convergence coefficients C for the energies of the low-lying 
bound states are shown. 



optimized for a given angular momentum. We report the estimates for the 
speed of convergence for s-, p-, d- and f-waves in Table [2j 

We have optimized the grid to reproduce 5 lowest eigenstates of the two- 
body Hamiltonian with angular momentum l opt and used those grids to cal- 
culate the ground state and low-lying vibrationally and rotationally excited 
states of the system. The extrapolated energy estimates (see the discussion 
for Eq. E]) agree excellently at the level of 10 -13 a.u., which is comparable 
with the error introduced by the cut-off distance of 500 a.u. As expected, 
optimizing the grid for the s-wave spectrum provides a very good basis for 
reproducing the rotationally excited states. Optimizing the grid for the ro- 
tationally excited states similarly gives very good results for the lowest two 
s-wave vibrational states. The convergence for the near-threshold s-wave 
state, however, is getting much slower, although still comparable with the 
speed of convergence for the low-lying states. We therefore conclude that only 
the s-wave optimization is needed to provide a good description of processes 
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involving the few lowest rotational excitations. 

The other important question is how many states should be included 
in the grading function fl5]). From Table [3] we can conclude that the very 
minimal number of states to be included in the grading function to reproduce 
the discrete spectrum of the operator should be equal to the number of the 
s-wave bound states. When including more states into the grading function 
the convergence of the bound state energies is getting insignificantly slower. 
We can expect, however, that the states from the continuum start being 
reproduced better. To check the convergence of the continuum states we fit 
the low-energy expansion for the s-wave scattering phase 

k cot 5 (k) = -i + ^+XA; 4 + 0(A; 6 ) (7) 

to the spectrum of the discretized Hamiltoniany and calculate the conver- 
gence coefficients of the effective range parameters, in quite the same man- 
ner as was done for the binding energy (Table H]). To estimate the scattering 
parameters we used the six lowest positive eigenvalues of the Hamiltonian. 
As there are three bound states in the system, it is not surprising that the 
fastest convergence is obtained with nine states included in the grading func- 
tion. Similar to the bound state case, including extra states into the grading 

2 Note that the energies of the discretized continuum are determined by the scattering 
phase and the boundary condition at the right end R of the interval. For the s-wave states 
the wave function behaves as sm(kR+S(k)), and the phase shifts can be recovered from the 
positive energy spectrum E n and the box size R using the condition \J E n R + 6(y/E n ) = 
nn. For large R the number of energies E n that fall in the region where the effective 
range expansion is valid becomes sufficient to extract the parameters of the effective range 
expansion. 
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Table 3: Ne2: speed of convergence for the s-wave vibrational bound states and first 3 
rotational excitations as a function of the number of optimized states. 



function does not have much effect on the speed of convergence after all the 
states used for estimating the effective range parameters have already been 
included. 

Before discussing the application to three-body calculations, we provide a 
check of the suitability of the grading function suggested in Ref . jlj] . For this 
purpose we compare the speed of convergence obtained with different values 
of k in Eq. but holding the order of the polynomial k! fixed, and check 
whether the value of k = k! provides the best speed of convergence. We sum- 
marize our observations in Tables |5l [6] and [7J In Table |5] we summarize our 
convergence speed analysis for the He2 dimer bound state energy (LM2M2 
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Table 4: Ne-Ne scattering parameters speed of convergence as a function of the number 
of states included in the grading function. Those data points marked with an asterisk 
indicate a failed convergence estimate. 
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potential |l4|) with the grid optimized for the ground and one continuum 
state. The theoretically optimal values are shown in bold font. In Table [6] 
we show the same data for the grids optimized for the bound state only. The 
lower values of k correspond to bigger variations of the mesh points density 
distributions, the bigger the k the more uniform the grid is. We can easily 
observe that for the values of k considerably different from the optimal val- 
ues the speed of convergence essentially deteriorates. For the values of /c a 
little bit smaller than the theoretically optimal one, however, the observed 
speed of convergence can be close to the optimal or even demonstrate a little 
better convergence. In the case of two states included into the grading func- 
tion (Table |5]) this can be partially attributed to the fact that the grading 
function is actually optimized for more than one state which we are track- 
ing. However, we observe similar behavior even if we optimize the grid to 
reproduce the bound state only. As our grids are constructed on the base 
of an iterative numerical procedure which is not proved to be exact and the 
speed of convergence parameter is intended as a qualitative measure of the 
convergence properties, this minor deviation of the observed optimal k from 
the theoretically optimal value can not be considered evident (see the discus- 
sion for Eq. ([()])). It is also useful to note that when the value of k exceeds 
the theoretically optimal value, the coefficient C starts to grow rapidly. 

In Table [7] we report the speeds of convergence for the energies of the s- 
wave bound states of the Ne 2 dimer. We have optimized the grid to reproduce 
the three bound states of the dimer and one state from the continuum. The 
speed of convergence coefficient for the near-threshold state clearly reaches 
its minimum at the theoretically optimal value of k, which confirms the 
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S'3,2 "$5,3 



k 


Ei (extrapolated) 


C 


E2 (extrapolated) 


C 


6 


-4.146595e-09 


2e-4 


n/a 


n/a 


7 


-4.146597e-09 


2e-5 


n/a 


n/a 


8 


-4.146598e-09 


-le-5 


n/a 


n/a 


9 


-4.146596e-09 


-5e-5 


-4.146638e-09 


-2e-2 


10 


-4.146598e-09 


-le-4 


-4.146578e-09 


3e-4 


11 


-4.146598e-09 


-3e-4 


-4.146577e-09 


8e-4 


12 


-4.146608e-09 


-6e-4 


-4.146579e-09 


le-4 


13 


-4.146636e-09 


-le-3 


-4.146581e-09 


6e-4 


14 


-4.146658e-09 


-le-3 


-4.146573e-09 


-5e-4 


15 


-4.146751e-09 


-9e-4 


-4.146586e-09 


2e-3 


16 


-4.146510e-09 


-le-2 


-4.146581e-09 


-7e-3 



Table 5: Optimality check of the automatically generated grids: He2 bound state energy 
(a.u.) and the speed of convergence for different values of k in Eq. ((5]). The grading 
function is optimized for n = 2 states. Bold font indicates the asymptotically optimal 
choice 
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53,2 




S$,3 




k 


Ei (extrapolated) 


C 


E 2 (extrapolated) 


C 


6 


-4.146595e-09 


2e-4 






7 


-4.146598e-09 


2e-5 






8 


-4.146595e-09 


-le-5 






9 


-4.146598e-09 


-6e-5 


-4.146589e-9 


-le-2 


10 


-4.146597e-09 


-2e-4 


-4.146622e-9 


8e-4 


11 


-4.146589e-09 


-3e-4 


-4.146634e-9 


8e-4 


12 


-4.146625e-09 


-7e-4 


-4.146624e-9 


4e-4 


13 


-4.146597e-09 


-2e-3 


-4.146631e-9 


4e-4 


14 


-4.146551e-09 


-2e-3 


-4.146642e-9 


-6e-4 


15 


-4.146658e-09 


-2e-3 


-4.146650e-9 


2e-3 


16 


-4.146452e-09 - 


-6e-3 


-4.146649e-9 


-3e-3 


17 


-4.146679e-09 


-le-2 


-4.146655e-9 


7e-3 


18 


-4.146670e-09 


-le-2 


-4.146656e-9 


-2e-2 


19 


-4.146366e-09 


-6e-3 


-4.146574e-9 


-7e-2 


20 


-4.146584e-09 


3e-2 


-4.146598e-9 


-3e-l 



Table 6: Optimality check of the automatically generated grids: He2 bound state energy 
(a.u.) and the speed of convergence for different values of k in Eq. {2}. The grading 
function is optimized for the bound state only. Bold font indicates the asymptotically 
optimal choice 
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5*3,2 






5*5,3 






k 




c Ne * 




CNe 2 






6 


-2e-5 


-2e-5 


3e-3 


n/a 


n/a 


n/a 


7 


-2e-4 


-5e-4 


5e-4 


n/a 


n/a 


n/a 


8 


-6e-4 


-le-3 


le-4 


n/a 


n/a 


n/a 


9 


-le-3 


-2e-3 


4e-6 


-3e-2 


8e-2 


-2e-2 


10 


-2e-3 


-3e-3 


-5e-5 


-3e-2 


le-2 


le-4 


11 


-4e-3 


-5e-3 


-8e-5 


-9e-2 


3e-3 


2e-3 


12 


-3e-3 


-5e-3 


-3e-5 


-2e-l 


-6e-3 


9e-4 


13 


2e-3 


2e-3 


2e-4 


-7e-l 


-3e-2 


-7e-4 


14 


2e-2 


2e-2 


9e-4 


-2 


-9e-2 


-4e-3 


15 


2e-2 


2e-2 


le-3 


-4 


-2e-l 


-le-2 


16 


5e-3 


4e-2 


2e-3 


-4 


-4e-l 


-2e-2 



Table 7: Optimality check of the automatically generated grids: the speed of convergence 
for the Ne2 dimer s-wave states. Bold font indicates the asymptotically optimal choice [lj]. 
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practicality of our approach. 



4. The three-body calculations with optimized grids 

As mentioned in the introduction, the practical value of optimizing the 
grids for two-body states lies in their application to more complex few-body 
calculations. In this case, reducing the number of points needed to achieve the 
required accuracy is critical for saving computational time. The direct use of 
the two-body optimized grids described in the previous section is especially 



appropriate in the Faddeev (three-body) [ll|, [12|, [13| or Faddeev-Yakubovsky 
(four-body) j^, O] formalism. When solving the Faddeev equations, as we 
shall briefly discuss below, it is natural to represent the grid as a direct prod- 
uct of the grid supporting the two-body bound states, and grids describing 
other coordinates in the configuration space of the three-body system. We 
shall demonstrate, that optimizing the grid even in one of the coordinates 
can improve the accuracy of the three-body calculations substantially. 

Here, for simplicity, we shall discuss only the three-body calculations 
below the three-body threshold, which physically means that we restrict our 
attention to three-body bound states or scattering with only two clusters in 
the initial and the final states of the system. Within the Faddeev approach 
the wave function of such states resolves into a sum of three components 

that correspond to different partitionings of the three-body system into an 
interacting two-body subsystem and the free third particle. Xi and yi are the 
Jacobi coordinates for the i-th partitioning: Xj connects the particles in the 
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i — th pair, y { points to the remaining particle from the 2-th pair center of 
mass. The asymptotic properties of the components <3>j below the three-body 
threshold are very simple: at large distances between the i-th particle and 
the center of mass of the corresponding two-body cluster {\yj\ — \ oo), the 
Faddeev components factorize. In the very simplest case of one single s-wave 
bound state in the pair the components <3>j behave as 

^ ~ 2 (W)/(N) , 

where 2 (|#i|) is the wave function of the two-body subsystem and f(\yi\) 
describes the free motion of the third particle. When solving the Faddeev 
equation numerically, we, therefore, need to be able to reproduce the behav- 
ior of the two-body clusters accurately. For this purpose we use the grid 
optimization procedure described above. In principle, it is also possible to 
develop some simple criteria for optimizing grids in the reaction coordinates 
y. This problem, however, is beyond the goals of this paper and we shall 
discuss it elsewhere. 

In order to demonstrate the importance of using optimized two-body 
grids in three-body calculations we have performed a series of calculations 
of He 3 ground and excited states with both optimized and non-optimized 
grids using a different number of points in the cluster coordinates. The 
number of points as well as the grid shape in the reaction coordinate y 
has been fixed for all sets of calculations. We used 100 grid points with 
S 5t3 splines in the y coordinate, the cutoff radius has been set to y ma x = 
2000 x \/3/2 a.u. and xitl mapping has been used to construct the non- 
uniform grid. As our three-body calculations here are only for purposes of 
demonstration, we used the simplest possible grid consisting of one single 
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interval in the angular coordinate (which is the cosine of the angle between 
the direction to the third particle and the axis of the He2 cluster). We used 
S§ t 3 splines in the angular coordinate. In this case the angular basis reduces 
to polynomials of 5-th degree. The angular symmetry of the system has 
not been taken into account explicitly, and, effectively, such angular basis 
corresponds to taking into account two first virtual rotational excitations 
(/ = 0, 2, 4). More details on the method we use to solve Faddeev equations 



can be found in 



11 



12 



131 ] . The results are shown in Figs. H] and EJ where 
we have plotted the energies of the trimer bound states calculated using S3 2 
(Fig. @J and 5*5,3 splines (Fig. [5]) as functions of the (appropriately scaled) 
number of grid points in the cluster coordinate. In both cases we observe 
much faster convergence and much smaller variation of the numerical results 
when using the optimized grids in cluster coordinates. We obtained the 
following estimates of the 4 He3 bound state energies. With cubic splines and 
100 grid points in x we find E = -3.98756 x 10" 7 a.u = -125.917 mK 
for the ground state and E x = -7.1984 x 10~ 9 a.u. = -2.2731 mK for 
the excited state. With quintic splines and 70 grid points the corresponding 
results are E = —3. 9877 x 10~ 7 a.u = —125.920 mK for the ground state and 
Ei = —7.1978 x 10 -9 a.u. = —2.2729 mK. These results agree perfectly well 
with previously reported independent results using the same angular basis 
but a different projection operator) Eq = —125.9 mK and E\ = —2.28 mK 
3). They also agree with previously reported results (for the same angular 



basis) of one of the authors 11] Eq = —125.81 mK and Ei = —2.2677 mK 
obtained with a similar method using less dense manually fine-tuned gridjfl 



A slightly different coupling constant employed in previous calculations [11] makes for 
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Figure 4: Energies of He3 bound states calculated using £3,2 splines as a function of the 
number of grid points in cluster coordinates 
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Figure 5: Energies of He3 bound states calculated using 53,2 splines as a function of the 
number of grid points in cluster coordinates. 
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5. Conclusions 



We have presented an approach to constructing an optimized nonuniform 
grid for use in quantum few-body calculations. The approach is based on 
the results of Ref. jl|, where a grading function which asymptotically mini- 
mizes the L2 norm of the interpolation error is introduced. We have slightly 
modified the optimization criterion to have several low-lying eigenstates of 
the two-body Hamiltonian interpolated well. 

We have studied the convergence properties of the optimized grids. For 
this purpose we have introduced the speed of convergence coefficient which 
characterizes the rate at which a physical observable - such as energy - 
converges as the number of grid points is increased. The optimized grids, 
even being only asymptotically optimal, demonstrate convergence properties 
superior to other choices of non-uniform grids routinely employed in few- 
body calculations. As far as rotational excitations of a two-body system can 
be taken into account perturbatively, it is sufficient to optimize the grids to 
reproduce the s-wave only and there is no need to optimize the grid for all the 
rotational excitations. The number of two-body states to be included into 
the grading function depends on the physical problem being solved. All the 
states that contribute to the long-range asymptote of the few-body system 
must be included into the grading function. Including some extra states can 
be beneficial, as it makes possible to account for the corresponding virtual 



60% of the discrepancy for the excited state and 100% of the discrepancy for the ground 
state value. The sensitivity of the results to the details of the interaction is a subject of a 
separate study being prepared for publication. 
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excitations more accurately, while the convergence of the low-lying states is 
only slightly affected. 

We have used the optimized grid in three-body calculations to describe 
cluster degrees of freedom. More accurate description of the internal dynam- 
ics of the colliding clusters leads to substantial improvements in the accuracy 
of calculations. The presented algorithm has allowed us to eliminate a dif- 
ficult and time-consuming stage of manual fine-tuning the nonuniform grids 
to be used with a particular three-body system. This result is essential for 
the goals of our project to develop a simple and effective public code for 
low-energy quantum three-body calculations. 

The presented algorithm has been constructed specifically for solving the 
quantum three-body problem with short-range interaction on the base of Fad- 
deev equations. The approach itself, however, has an unexplored potential. 
Let us make a few remarks on other possible applications. Optimal grids for 
Coulomb systems can be considered. In this case the grading function can 
be constructed analytically. Adiabatic hypersphecical calculations can also 
benefit from adding a grid adaptation step after calculating the adiabatic 
basis for each value of the hyperradius. Finally, finding an effective method 
to optimize the "reaction" degrees of freedom in Faddeev calculations is a 
subject especially interesting for our project. 
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